%% This script uses the CHAMP EEF and ACE IEF to find the transfer function between them

clear;

%% JULIA data
% load /data/backup/mnair/longp/alldays JULI_SEG ACE_SEG;
% JULI_SEG = JULI_SEG.*24.366*1e-3; %mV/m
% [Txy_short,F_short] = tfestimate(ACE_SEG,JULI_SEG,hanning(72),0,72,1/(5*60)); %Txy is agains calculated just to get
% [Cxy_short,F_short] = mscohere(ACE_SEG,JULI_SEG,hanning(72),0,72,1/(5*60));
% N_data = length(ACE_SEG) / 72;
% Err_short = sqrt( 1/(2*(N_data-1)) .* ( (1-Cxy_short)./Cxy_short ) ) .* abs(Txy_short);
% Err_short_lg = 0.434 * Err_short ./ abs (Txy_short);
%load /data/backup/mnair/longp/eef_data_mod.mat eef;

load /nfs/satmag_work/mnair/projects/longp/eef_data_2000_2009.mat eef;



% solar activity
% L = eef(:,1) <= datenum(2005,08,1);
% 
% eef(L,:) = [];

% 

% # Field 1: timestamp (UT)

% # Field 2: longitude (degrees)

% # Field 3: latitude (degrees)

% # Field 4: local time (hours)

% # Field 5: season (day of year)

% # Field 6: eastward electric field (mV/m)

% # Field 7: equatorial vertical electric field at 105 km altitude (mV/m)

% # Field 8: equatorial (UxB)_z 105 km altitude (mV/m)

% # Field 9: DC current shift (A/m)

% # Field 10: CHAMP peak current value (A/m)

% # Field 11: KP

% # Field 12: F10.7 (W/m^2)

% # Field 13: F10.7A (W/m^2)

% # Field 14: R^2 (coefficient of determination)

% # Field 15: R (correlation of CHAMP and model profiles)

% # Field 16: chi^2

% # Field 17: EEFM modeled eastward electric field (mV/m)



% I tried to see weather the TF magnitude plot can be reproduced by just

% two random variables.

% std_eef = std(eef(:,6));

% mean_eef = mean(eef(:,6));

% eef_rand = mean_eef + randn(length(eef),1)*std_eef;

% eef(:,6) = eef_rand;

% eef(:,17) = mean_eef + randn(length(eef),1)*std_eef;

%





%load /data/backup/mnair/longp/OMNI_ELEC_new ace_all;

load /nfs/satmag_work/mnair/projects/ace_tensor/acedata/ace_2000_2010.mat ace_all;

%ace_all(:,1) = ace_all(:,1) - datenum(2000,1,1); % Because the eef data uses fday wrt 2000



% Replace ace data with random signals



% L = isnan(ace_all(:,2)) | isnan(ace_all(:,3));

% ace_all(L,2) = 0;

% std_ace = std(ace_all(:,2));

% mean_ace = mean(ace_all(:,2));

% ace_rand = mean_ace + randn(length(ace_all),1)*std_ace;

% ace_all(:,2) = ace_rand;

% ace_all(:,3) = mean_ace + randn(length(ace_all),1)*std_ace;

%

% ace_all(L,2) = NaN;

% ace_all(L,3) = NaN;

%

%



load /nfs/satmag_work/mnair/projects/longp/aplist.mat;



% Initializations



options.ap_lower_limit = 20; % FIXED at 20 lower limt of Ap

options.ap_upper_limit = 1000;

options.phase_delay = 17;% FIXED at 17 time delay for IEF to propagate from Bow-shock to ionopshere minutes

options.plc = 'c';% plot color

%options.tol = 1.8/24; % max time interval tolerated

options.tol = 0.0750;

%options.min_time_length = 128/24; % minimum time length required in decimal days

options.min_time_length = (128*8)/24;

options.lt_start = 9; % FIXED 9-15 (data permits 7-17)

options.lt_end = 15;

options.overlap = 1;

options.champ_resampling = 1;

options.coh_plot_lim = 0.0;

options.ace_interp_method = 'spline'; %FIXED spline



options.minimum_period_plot = 6 ; % 5 if the min length is around 30-40 hours

options.maximum_period_plot = 512;



% options.minimum_period_plot = 18; % 5 if the min length is around 30-40 hours

% options.maximum_period_plot = 50;





options.Eeff_cut_off = 100; % The cut off amplitude for IEF .LUHR & MAUS EPS 2010 used a value of 8



options.des_int = 2; %desired sampling interval in hours

% the sampling interval is set ~30 minutes above the

% champ EQ revisit time (92 minutes)



options.ace_gap_tol = ( options.min_time_length + 6/24 ) * 0.1 ;

options.ace_gap_len_tol = floor ( options.ace_gap_tol * 1440/5 );

% the gap tolerence is 10% of the ace data length



% finding gaps in ACE data with gaps > 12 minutes



nd = 1;

np = 0;



L = isnan(ace_all(:,2)) | isnan(ace_all(:,3));



for i = 1: length(ace_all) - 1,

    

    

    if  L(i) && L(i + 1)

        

        np = np + 1;

        

    else

        if ( ace_all(i,1) - ace_all(i - np ,1) > 12/1440 )

            % Mar 2012. One of the issues to be looked at is how these

            % gaps are affecting the spectra. A future improvement

            % should be to break the ACE data into pieces (just like

            % EEF data) and limit the Fourier analysis to within the

            % segements.

            data_index_ace(nd,2) = i;

            data_index_ace(nd,1) = i - np;

            nd = nd + 1;

            

        end;

        np = 0;

        

    end;

end;





%%



% Calculating the effective electric field (LUHR EPS)

% E = 8IEF sqrt (64 + IEF^2). This has the effect of restructing amplitudes

% > Eeff_cut_off to within Eeff_cut_off. The reason to introduce this was

% a small number of high energey events (< 2% of the data) can cause the issues with th eLS

% welch power spectra estimation. I haven't looked at how this is going to

% affect the powe spectra though. - When I looked the ACE power spectra,

% this amplitude modulation as an effect of reducing the power (more at

% longer periods). This will give a false raise in TF amplitude. So better

% not use it. ( I have done this by increasing the Eeff_cut_off to 100 )



% The natural power spectra of IEF (and most geomagnetic variations) obey

% the slope -2 rule. That is gain of two orders of magnitude for every

% order of frequency decrease. Hence, the larger amplitude variations on

% average and in general correspond to the long period fluctuations.In

% other words, the larger the amplitude fluctuations, the larger the

% possibility that it is a long period one. Hence by selectively depressing

% the large amplitude signals, we disturb the slope -2 rule. Hence do not

% use the following scripts for E-eff 



% I hope this what it was happening on the IEF power spectrum. 

L = ace_all(:,2) < 0;

temp = abs(ace_all(:,2));

Eeff = options.Eeff_cut_off * temp ./ sqrt (options.Eeff_cut_off^2 + temp.^2);

ace_all(:,2) = Eeff;

ace_all(L,2) = ace_all(L,2) * -1;

% 

%IEF_Ez



L = ace_all(:,3) < 0;

temp = abs(ace_all(:,3));

Eeff = options.Eeff_cut_off * temp ./ sqrt (options.Eeff_cut_off^2 + temp.^2);

ace_all(:,3) = Eeff;

ace_all(L,3) = ace_all(L,3) * -1;





% down sampling and resampling to original interval



%The following script is to resample the ace data at a lower rate.

% The matlab resample apply a low-pass filter before resampling

% this allows for some remedy to the aliasing issues. This is a very

% important step. With out this sever aliasing issues can creep into the

% spectra analyssis.



L = isnan(ace_all(:,2));

y = interp1(ace_all(~L,1),ace_all(~L,2),ace_all(:,1),options.ace_interp_method);

ace_down = resample(y,1,options.des_int*60/5); % Checked the resampling time axis OK

ace_inter = interp1(ace_all(1:options.des_int*60/5:end,1),ace_down,ace_all(:,1), options.ace_interp_method);

ace_all(:,2) = ace_inter;





%IEF Ez



L = isnan(ace_all(:,3));

y = interp1(ace_all(~L,1),ace_all(~L,3),ace_all(:,1),options.ace_interp_method);

ace_down = resample(y,1,options.des_int*60/5); % Checked the resampling time axis OK

ace_inter = interp1(ace_all(1:options.des_int*60/5:end,1),ace_down,ace_all(:,1), options.ace_interp_method);

ace_all(:,3) = ace_inter;



% make sure that the ACE data gaps > 12 min are filled with zeros using the

%% data index determined above





for i = 1:length(data_index_ace),

    

    %     ace_all(data_index_ace(i,1):data_index_ace(i,2),2) = 0;

    %     ace_all(data_index_ace(i,1):data_index_ace(i,2),3) = 0;

    %

    %linear interpolation for the gaps > 12 minutes. Note that for gaps

    %less than 12 minutes, we use the defined interpolation method (such as

    %spline)

    ace_all(data_index_ace(i,1):data_index_ace(i,2),2) = ...
        interp1(ace_all([data_index_ace(i,1)-1 data_index_ace(i,2)+1],1), ...
        ace_all([data_index_ace(i,1)-1 data_index_ace(i,2)+1],2), ...
        ace_all(data_index_ace(i,1):data_index_ace(i,2),1));

    

    ace_all(data_index_ace(i,1):data_index_ace(i,2),3) = ...
        interp1(ace_all([data_index_ace(i,1)-1 data_index_ace(i,2)+1],1), ...        ace_all([data_index_ace(i,1)-1 data_index_ace(i,2)+1],3), ...
ace_all(data_index_ace(i,1):data_index_ace(i,2),1));

    

    

    

end;





options.ace_gap_tol = ( options.min_time_length + 6/24 ) * 0.1 ;

options.len = floor( options.min_time_length *24 / options.des_int ); % This number is the time length (also data length) in hours

lt_incr = 1;          % An index used for counting the LT window numbers







%ace_all= [ace_all(1:des_int*60/5:end,1) ace_down];



% Arrange the CHAMP data into segments

% the climatology is limted to > 7 and < 17 LT


%% Select CHAMP data


nd = 1;

np = 0;

% options.lt_start = 7;
% options.lt_end = 17;

for i = 1: length(eef) -1,

    

    if (eef(i+1,1) - eef(i,1) <= options.tol && eef(i,4) > options.lt_start && eef(i,4) < options.lt_end )

%             if (eef(i+1,1) - eef(i,1) <= options.tol && eef(i,4) > options.lt_start && eef(i,4) < options.lt_end ...
%             && eef(i,1) - eef(i - np ,1) <= options.min_time_length )


        np = np + 1;

    else

        % if (np >=16)

        if ( eef(i,1) - eef(i - np ,1) >= options.min_time_length )

            data_index(nd,2) = i;

            data_index(nd,1) = i - np;

            np = 0;

            nd = nd + 1;

        else

            np=0;

            

        end;

    end;

    

end;





%%

% If overlap is enabled, change the data_index to ensure a

% 50 % overlap between the windows. This should be done only for continous

% windows. The overlapped windows are added to the end of the data_index. 



if options.overlap == 1

    

    nd = 1;

    np = 0;

    lendata = length(data_index) ;

    

    for i = 1:lendata - 1

        

        if data_index(i+1,1) - data_index(i,2) == 1,

            data_index(lendata + nd, 1) = floor(sum(data_index(i,:))/2);

            data_index(lendata + nd, 2) = floor(sum(data_index(i+1,:))/2);

            nd = nd + 1;

        end

        

    end;

end;



% iterations

CHAMP_SEG = [];

ACE_SEG = [];

TIME_SEG = [];



N_seg = 1;%the increasing counter for array JULI_SEG & ACE_SEG

N_data = 0;



for i = 1: size(data_index,1),

    

    

    L = ace_all(:,1) + options.phase_delay/(60*24) >= eef(data_index(i,1),1) - 3/24 ...
        & ace_all(:,1) + options.phase_delay/(60*24) <= eef(data_index(i,2),1)+ 3/24 ;

    

    if sum(L) > 0,

        ace_time = ace_all(L,1) +  options.phase_delay/(60*24) ;

        ace_data = ace_all(L,2);

        ace_data_ez = ace_all(L,3);

        L = isnan(ace_data);

        if any(ace_data),

            if  max(diff([ace_time(1) ; ace_time(~L) ; ace_time(end)] )) < options.ace_gap_tol && ...
                    sum(L) < options.ace_gap_len_tol,

                

                



                if abs(mean(diff(ace_time))-0.0035) <= 1e-004, %Use this with ACE 5 min averages



                    

                    dummy = eef(data_index(i,1):data_index(i,2), 6 );% - eef(data_index(i,1):data_index(i,2), 17 );

                    time_axis_desired = eef(data_index(i,1), 1 ):options.des_int/24:eef(data_index(i,2), 1 );

                   

                    

                    %                   Resampling the CHAMP data from the average sampling

                    %                   interval (about 92 minute) to 120 minutes using the

                    %                   resample function.

                    % get a rational number for the original/desired

                    % sampling ratio

                        if options.champ_resampling == 1,

                                        [n,d] = rat(median(diff(eef(data_index(i,1):data_index(i,2),1)))/(options.des_int/24));

                                        y = resample(dummy,  n,  d);  % Now resample it

                                        t2 = (0:(length(y)-1)) * d / (n * 1/ median(diff(eef(data_index(i,1):data_index(i,2),1))) );

                                        CHAMP_EEF = interp1(t2 + eef(data_index(i,1),1), y, time_axis_desired);

                        

                        else,

                        

                             CHAMP_EEF = interp1(eef(data_index(i,1):data_index(i,2),1), dummy,time_axis_desired);

                        

                        end;

                    

                    

                    

                    ACED = interp1(ace_time, ace_data, time_axis_desired);

                    ACEZ = interp1(ace_time, ace_data_ez, time_axis_desired);

                    L = isnan(ACED);

                    ACED(L) = [];

                    

                    L = fday_ap >= eef(data_index(i,1), 1 )& fday_ap <= eef(data_index(i,2), 1 );

                    mean_ap = mean(ap(L));

                    

                    if sum(isnan(ACED)) < 1 && mean_ap >= options.ap_lower_limit ...
                            && length(ACED) >= options.len && length(CHAMP_EEF) ... 
                            >= options.len && mean_ap <= options.ap_upper_limit;

                                                

                        CHAMP_SEG(N_seg:N_seg+(options.len-1)) = CHAMP_EEF(1:options.len);

                        ACE_SEG(N_seg:N_seg+(options.len-1)) = ACED(1:options.len);

                        ACE_SEG_EZ(N_seg:N_seg+(options.len-1)) = ACEZ(1:options.len);

                        TIME_SEG(N_seg:N_seg+(options.len-1)) = time_axis_desired(1:options.len);

                        % IMF_BZ_SEG(N_seg:N_seg+71) = IMFBZ(1:72);

                        N_seg = N_seg+options.len;

                        N_data = N_data+1;

                        

                        

                    end;

                   



                else,

                    %fprintf('Day %d has some missing time stamp\n', Julia_W(i).fday);

                end;

            end;

        end;

    end;

end;



fprintf('Done ! N_data = %3d , Length = %d hours, LT = %d-%d, AP limit = %d color = %s, interp = %s\n', ...
    N_data, options.min_time_length*24, options.lt_start,options.lt_end, options.ap_lower_limit, options.plc, options.ace_interp_method);



%%

% Coherence, phase and tranfer function

%%JULI_SEG = JULI_SEG.*24.366*1e-3; %mV/m



open('/nfs/satmag_work/mnair/projects/longp/Figures/julia_with_error_bar.fig');

hold on

[Cxy_long,F_long] = mscohere(ACE_SEG, CHAMP_SEG,hanning(options.len),0,options.len,1/(options.des_int*3600)); %1/(5*60) = sampling frequency in Hz )

[Pxy,F_long] = cpsd(CHAMP_SEG,ACE_SEG,hanning(options.len),0,options.len,1/(options.des_int*3600)); %

[Pxx,F_long] = pwelch(CHAMP_SEG,hanning(options.len),0,options.len,1/(options.des_int*3600));

[Pyy,F_long] = pwelch(ACE_SEG,hanning(options.len),0,options.len,1/(options.des_int*3600));

phase = angle(Pxy);

[Txy_long,F_long] = tfestimate(ACE_SEG,CHAMP_SEG, hanning(options.len),0,options.len,1/(options.des_int*3600));

%tf = conj(Txy');

L = Cxy_long > options.coh_plot_lim & (1./F_long)./(3600) >= options.minimum_period_plot &...
    (1./F_long)./(3600) <= options.maximum_period_plot ;



Err_long = sqrt( 1/(2*(N_data-1)) .* ( (1-Cxy_long)./Cxy_long ) ) .* abs(Txy_long);



Err_long_lg = 0.434 * Err_long ./ abs (Txy_long);



errorbar(log10((1./F_long(L))./(3600)), log10(abs(Txy_long(L))),(Err_long_lg(L)),(Err_long_lg(L)),options.plc, 'linewidth',2);

hold on;

axis([ -inf   inf   -3   -1]);



for i = 2:length(Cxy_long),

    fprintf('%06.3f ',(1./F_long(i))./(3600));

end;

fprintf('\n');

for i = 2:length(Cxy_long),

    fprintf('%06.3f ', Cxy_long(i));

end;

fprintf('\n');

%%

% title_str = sprintf(' Proc%3d-%dLT%d-%dAP%d%s\n',N_data, options.min_time_length*24, options.lt_start,options.lt_end, options.ap_lower_limit, options.ace_interp_method);

% title(title_str);



% The CHAMP EEF spectra has peaks at 6, 8, 12 and 24. While this is not

% surprising, the 6 hour harminics has higher amplitudes than 8 and 12.

% Checking this with fitting sinusoids to the data



    

% periods = [4 4.8 6 8 11.967236 12  12.421 12.6583 23.934472 24 25.891];

% periods = sort(unique([1./(3600*F_long(2:end))' 4 6 8 12 24])) ;

% t = TIME_SEG;

% 

% x = [cos(2*pi*t/( periods(1)/24))'       sin(2*pi*t/( periods(1)/24))'];

% 

% for j = 2:length(periods),

%     

%     x = [x cos(2*pi*t/( periods(j)/24))'       sin(2*pi*t/( periods(j)/24))'];

%     

% end;

%     

% spectra_rob_eef = robustfit(x,CHAMP_SEG);

% spectra_data_rob_eef = abs(complex(spectra_rob_eef(2:2:end),spectra_rob_eef(3:2:end)));

